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Recent investigations have provided important insights into the complex structure and dynam¬ 
ics of collectively moving flocks of living organisms. Two intriguing observations are - scale-free 
correlations in the velocity fluctuations, in the presence of a high degree of order, and topological 
distance mediated interactions. Understanding these features, especially, the origin of fluctuations, 
appears to be challenging in the current scheme of models. It has been argued that flocks are 
poised at criticality. We present a self-propelled particle model where neighbourhoods and forces 
are defined through topology based rules. The force fluctuations occur spontaneously, and gives 
rise to scale-free correlations in the absence of noise and in the presence of alignment of velocities. 
We characterize the behaviour of the model through power spectral densities and the Lyapunov 
spectrum. Our investigations suggest self-organized criticality as a probable route to the existence 
of criticality in flocks. 

PACS numbers: 87.18.-h,5.65.-|-b,5.50.+q,75.10.Hk 


The last decade has seen a growing interest in the study 
of collective motion of living as well as non-living sys¬ 
tems. Examples include, collectively swimming bacteria 
[1], human moshing [2], propulsion of Janus particles [3] 
and motion in granular systems [4]. Particularly worth 
mentioning are the investigations of the StarFlag project 
on the structure and dynamics of three dimensional star¬ 
ling flocks. It was found that any bird in a flock of star¬ 
lings, on the average, interacts with a fixed number of its 
nearest neighbours [5]. This notion of ‘topological dis¬ 
tance’, for defining the neighbourhood, is contrary to the 
assumptions in models [6] that a bird would interact with 
all other birds within a fixed metric distance. Moreover, 
the analyses revealed the existence of scale-free correla¬ 
tions in the velocity fluctuations of the birds [7]. Such 
scale-free correlations were also found in motile bacterial 
colonies [8] and midge swarms [9]. Interestingly, in case 
of starlings and bacteria the flocks showed a degree of 
order. The velocities of the individuals were found to be 
aligned, on the average. 

Critical fluctuations in models are generally obtained 
by tuning control parameters like noise and fields. How¬ 
ever, the mechanism through which the fine tuning of the 
parameters transpires in real flocks is not yet known. It 
has been argued that flocks are actually poised at critical¬ 
ity [10]. Maximum entropy modeling based on data [II], 
has been used to extract the values of the parameters that 
make a flock critical. The contrasting paradigm of self- 
organized criticality (SOC) has also been suggested as a 
dynamical route to the realization of criticality [7, 10]. 
In this paper, we provide a self-propelled particle (SPP) 
model along the latter lines. In our model, the critical 
state is realized in the absence of any external noise or 
field. We implement a topological distance based rule 
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for the choice of the neighbourhood of the particles. We 
assume that particles in the flock interact with a cer¬ 
tain fixed number nearest-neighbours through spring-like 
forces. We assume that particles have a preference to 
get attracted towards the neighbours that are located to¬ 
wards the centre of mass (COM). Previously, one of the 
authors considered a comoving boundary, defined with 
respect to the COM, where, the particles at or outside 
the boundary would be attracted towards the COM [12]. 
In a recent investigation, employing Vornoi tessellation, 
the particles at the convex hulls of the surface were con¬ 
sidered to get attracted towards bulk [13]. 

We consider N identical particles on the two dimen¬ 
sional space characterized by the positions, {rj}, and ve¬ 
locities {iTi}: where, 1 < i < N. We assume that the 
particles interact with each other through forces depend¬ 
ing on topological distances. At every instant, the mo¬ 
tion of a particle is influenced by its 2m “neighbours”, 
where m is a parameter in the model. To decide the in¬ 
stantaneous set of neighbours for particle i we divide the 
entire space into two non-overlapping regions, using the 
position of i and the location of the COM for the entire 
set of N particles. We consider the particles lying in the 
region, in which the COM lies, as particles located in 
the direction of the COM, with respect to the particle 
i. The m nearest particles in this region constitute the 
set Qi- Similarly, the m nearest particles in the other 
region make up the set Pi- The sets. Pi and Qi, together 
comprise the neighbourhood of i. This rule is illustrated 
in Fig. I. 

At every instant of time the sets Q and P for each 
of the particles are constructed and the positions and 
velocities are calculated using the equations of motion: 
dvijdt = Vi, and, M dvi/dt = F{vi) + FQ,i + Fp,i, where, 
M is the mass of each particle. The propulsion force, 
F{vi) = —b(^\vf\/vQ — l)vi, is the Rayleigh-Helmholtz 
(RH) friction [14-16], where, vq is the optimal speed. 
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FIG. 1. (a) The figure illustrates the choice of neighbours for 
a particle i with m = 3 at a certain instant of time. The 
dotted line connects the location of i to the centre of mass 
of the system (COM, denoted by a cross). The dashed line 
Li is perpendicular to the dotted line and passes through 
the location of i. The line Li divides the two dimensional 
space into two non-overlapping regions. We first consider the 
region in which the COM is located. We consider the particles 
lying in this region as particles located in the direction of the 
COM. (In general, particles lying on Li are also included in 
this region). Such particles are shown as filled black circles. 
In this region, the particles qi .2 and gi ,3 are the first 
three (m = 3) nearest neighbours of i and constitute the 
set Qi. The particles in the other region (denoted by filled 
grey circles) are considered as particles located in the direction 
opposite to the direction of the COM. The set of neighbours 
in this region is Pi — {pi,i,Pi, 2 ,Pi, 3 }. Together the sets Qi 
and Pi constitute the neighbourhood of i and are updated at 
each instant. The figures (b) and (c) schematically shows the 
sets Q and P for individual particles when m = 1 in a system 
of 20 particles at a certain instant. In the figure (b) an arrow 
points from each particle to its neighbours in the set Q. In (c) 
the arrow connects to neighbours in P. Note, that for some 
of the particles there are no neighbours in the set P. 


which makes the propulsion force zero; and the parame¬ 
ter b controls the strength of the propulsion force relative 
to the other forces. The term Fq^i and Fp^i are the forces 
acting on i due to the neighbours in the sets Qi and Pi, 
respectively. We assume that the force of interaction of i 
with a neighbour j in the set Qi is essentially spring-like 
with pQ^i = -{l/riQ^i) J2jGQi k{\rij\- a)rij/\rij\, where, 
Pj = Ti — Tj, k is the spring constant, a is the optimal 
separation with a neighbour, and ng^i is the number of 
neighbours in the set Qi. However, for the neighbours in 
Pi we consider only repulsive interactions below the sep¬ 
aration a, such that, Fp^i = —(1/np^i) 

I 

a)rij/\rij\, where, fij = I when |r) - r,| < a and /„ = 0, 
otherwise. We note here that, the condition ng^i = m is 
always satisfied when m N, but this is not the case 
with np^i. These can be seen in Fig. 1(b) and (c). For 
such particles the total number of neighbours is less than 
2m. 

We study the time evolution of the model by nu¬ 
merically integrating the equations of motion using the 
velocity-Verlet algorithm [17]. We set the scales for mea¬ 
surement as M, a and vq. All quantities are measured in 


these scales. For our main investigation, we take 6 = 1.0 
and k = 3.0, and the integration time step dt = 0.005. 
At every time step we update the set of neighbours for 
every particle. However, in the absence of a cut-off dis¬ 
tance for the cohesive force, the construction of sets P 
and Q becomes computationally intensive. 

At the beginning of the simulations we randomly po¬ 
sition all the particles in circular a region of having ra¬ 
dius '/N /2 and with randomly oriented velocities having 
magnitude unity. We simulate the model in open bound¬ 
ary conditions and, in general, we find that the flock is 
a single cohesive entity. On the average the particles re¬ 
main bounded in a finite region when in the frame of the 
COM. This is ensured by cohesive forces on the particles 
in the direction of the COM as well as the absence of any 
external noise. Although, artificial initial conditions for 
positions and velocities can be constructed which may 
make the flock to expand in space as time progresses, 
we did not encounter such cases, starting from arbitrary 
initial conditions. 

The velocity of the COM for the N particles is, V = 
Tf Si P- We characterize the alignment of velocities us¬ 
ing the speed, V = \V\. We would like to mention 
that, having attractive interactions with neighbours in 
P, causes the flock to stretch and with H ~ 0 producing 
weak translation of the COM. We exclude such interac¬ 
tions in our model. In Fig. 2(a) we show the approach 
to the steady state in our model, for five different values 
m. The values of V in the steady state is plotted in Fig. 
2(b). In general, for values of m > 3 the model produces 
flocks having V ~ 1.0. The case m = 3 is shown in Fig. 
3(a). For m < 3, the flocks fail to show any global order. 
Remarkably, in the case m = 3, not only the value of 

V is high, the fluctuations in V are high as well (from 
Fig. 2 (b), V = 0.977 ± 0.025). This also results in 
the tortuous trajectory for the COM (Fig. 2(c)). This 
is in contrast with the cases, m = 4 and 5, where the 
fluctuations appear to be rather suppressed (for m = 5, 

V = 0.997 ± 0.002). This is the main motivation for 
choosing m = 3 for the investigation of correlations in 
velocity fluctuations. 

The fluctuation in the velocity of a particle i, is defined 
as, Ui = Vi — V. In Fig. 3(b) we plot the fluctuations in 
the velocities for the particles when m = 3. We define 
the linear size (A) of the flock in the following fashion. At 
any instant we find the size by calculating the maximum 
of all distances jr) —rjj, where, i,j G [1, A^] [7, 8]. We find 
L by averaging the size over initial conditions and time 
snapshots. In Fig. 2(d) we plot L for different values of 
m. The plot reveals larger variation in the shape of the 
flock for smaller values of m. For, m = 3, we consider 
the variation in the size of the flock with the number of 
particles in system. In Fig. 2(e) we plot L versus W. We 
find that L ~ with = 0.5. For bacterial clusters 
moving in two dimensions this exponent was found to be 
0.6 [8]. Starling flocks have three-dimensional structure 
and have been found to be asymmetrical. Interestingly, 
a regression fit, to the data provided [7] for the number 
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FIG. 2. (Color online) (a) The time evolution of the configu¬ 
ration averaged speed of the centre of mass (V) is shown for 
different values of m. (b) The value of V in the steady state 
(as seen from the figure (a)) is plotted against m. (c) A typ¬ 
ical trajectory of the centre of mass for a flock with m — 3, 
monitored over 2000 time steps, (d) The linear size of the 
flock (L) in the steady state is plotted against m. For (a)- 
(d), N = 512. The straight lines in (b) and (d) are guides to 
the eye. (e) Log-log plot of the linear size of the flock (L) plot¬ 
ted against number of particles {N) when m = 3. The dashed 
line is suggestive of L increasing as with ("i = 0.5. The 
bars indicate standard deviation. 



FIG. 3. (a) A typical flock in the steady state generated 

by the model corresponding to N — 200 and m = 3. The 
instantaneous velocity vectors of the particles are plotted, (b) 
The fluctuations in velocities (with respect to the velocity of 
the COM) corresponding to the figure (a) are plotted. 


of birds in flocks and the corresponding sizes, results in 
the exponent 0.6. 

The spatial correlation function for veloc¬ 
ity fluctuations [7, 18] is defined as, C{r) = 

(l/C'o)Eij where, 

Co is a normalization constant, such that, C'(O) = 1. 
The Dirac ^-function is realized through the binning 
of the correlation function. The correlation length, 
is defined through the relation C{r = ^) = 0, and is 
expected to provide the average size of the correlated 
domains. For starling flocks [7] and bacterial clusters [8], 
the correlation length is found to be proportional to the 
linear size (L) of the system. Absence of a characteristic 


scale for domains implies scale-free correlations. A 
finite-size scaling ansatz for the correlation function 
shows that the correlation function evaluated in the limit 
of T —?► oo, C'oo('r) ^ r~'^. The exponent is estimated 
from the relation C'{r/^ = 1) ~ —L~^ [7]. Starling 
flocks revealed a extremely long-ranged correlations 
(7 ~ 0). Models have shown the existence of scale-free 
correlations in the highly ordered regime [18-20] mainly 
through the relation, ^(T) oc L. A 2d SPP model [21] 
using topology (Voronoi cells) to determine nearest 
neighbours, studied in the low-noise regime, yielded 
7 = 0.4 [18]. The 3d classical Heisenberg model with a 
dynamical external field and at low noise, showed 7 is 1 
and 0, in weak and strong fields, respectively [18]. 

We plot C'(r) in Fig. 4(a) for nine different values of 
N. We calculate ^ by averaging the positions of zero- 
crossings of C{r). The Fig. 4(b) suggests that the cor¬ 
relations are indeed scale-free. We find that, the gen¬ 
eral behaviour for large values of L is, ^(L) Ri ci.L, 
with, Cl = 0.29. The values of Ci in the cases of star¬ 
ling flocks and bacterial clusters are 0.35 and 0.3, re¬ 
spectively. In Fig. 4(c) we plot the correlation func¬ 
tions against rescaled length r/^. The hgure suggests a 
decay of the slope calculated a r/^ = 1 with increase 
in N. We determined the slopes, lC"(r/^ = l)j during 
the averaging of C{r) and The absolute values of the 
slopes are plotted with N in Fig. 4(d). A power-law fit 
suggests lC'(r/^ = l)j ^ with (2 = 0.22. Com¬ 

bining this result with dependence of L on N, we find, 
7 = C2/C1 = 0.44. This value suggests that in our model 
the decay of correlations is faster compared to starling 
flocks but is very similar to that observed in the low- 
noise regime of the topological SPP model [18, 21]. We 
would like to mention here that we also investigated the 
correlations in the fluctuation of speed [7, 8]. We find 
the overall behaviour to be qualitatively similar to that 
of correlations in velocity fluctuations. However, the ab¬ 
solute values for correlations and correlation lengths ap¬ 
pear to be lesser. A possible reason for this would be the 
RH friction which suppresses the fluctuations in speed. 

We also study the diffusion of particles in the COM 
frame. For starling flocks the diffusive behaviour reflects 
the fact that the set of neighbours for a given bird keeps 
changing in time [22]. At any given instant we consider 
a particle to be at the boundary if the number of neigh¬ 
bours in the set P is less than m. We tag individual 
particles and follow their motion until they reach the 
boundary. We calculate the average of the squared dis¬ 
placements, ((Ar')^), of these particles in the centre of 
mass frame. The Fig. 5(a) reveals, ((Ar'(t))^)^ , with 

/3 = 1.09. This indicates a normal diffusive behaviour in 
contrast to super-diffusive behaviour observed for star¬ 
lings, where /3 = 1.73. 

We investigated the nature of the spontaneous fluctu¬ 
ations by analyzing the time variation of different quan¬ 
tities. We calculated the power spectral density (PSD) 
corresponding to the time series’ for the a:-component of 
the velocity {vx) of a particle and the cc-component of its 
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FIG. 4. (Color online) (a) Plot of the correlation functions 
C(r) for velocity fluctuations for flocks with different N. (b) 
Plot of correlation length ^ versus the linear size L. The 
dashed line is a fit of the form ^ = cq + ciL, giving cq = 1.17 
and Cl = 0.29. (c) C{r) plotted against rescaled separation 
r/^. The symbols denoting different values of N in (a) and 
(c) are provided in the legend of (a), (d) Log-log plot of the 
slopes of the correlation functions plotted in (c) at r/^ = 1, 
as a function of N. The dashed line is a power-law fit having 
the form |C'(r/5 = 1)| ~ N~‘^^ giving (2 = 0.22. The bars 
in (b) and (d) indicate the standard deviations obtained from 
averaging over 100 time snapshots (separated by 10"* time 
steps) in the steady state starting from around 100 different 
initial conditions. 


position {x') in the COM frame. The results are shown in 
Fig. 5(b). The PSDs are averaged over several time series 
and plotted against frequency uj. We find the PSD for Vx 
has a power-law region of the form , with, a^, = 1.08. 

The upper bound to the power-law approximately coin¬ 
cides with the natural frequency, corresponding to the 
spring-like interactions, that is y^k/M = -x/S (in units of 
Vo/a). The peak appearing in the lower frequency region 
can be be attributed to the direction reversals of the par¬ 
ticle near the boundary of the flock. The PSD for x' is a 
power-law over most of the frequency span and decays as 
with, ax = 2.21. This signifies Brownian nature of 
the time series for x' and confirms the diffusive behaviour 
seen from Fig. 5(a). We believe that the origin of scale- 
free correlations in the absence of noise suggests that the 
system is in a state of self-organized criticality (SOC). 
The l/w"” decay of the PSD for velocity reaffirms this 
idea. 

The reason that the fluctuations persist in the system 
when m > 3, even in the absence of noise, is mainly 
the absence of any stable structure in the COM frame. 
Consider the subsystem shown in Fig. 5(c) where to = 3. 
The figure shows a particle A and its neighbours in the set 
Qa- We assume that any at any instant Qa = {B, C, D} 
and the distances between A and the neighbours is just 
greater than a, the optimal separation. Also, let all the 
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FIG. 5. (a) Log-log plot of the configuration averaged mean 
square displacement in the COM frame for tagged particles 
against time. The different symbols correspond to different 
flock sizes as indicated in the legend. The straight line, hav¬ 
ing slope /3 = 1.09, is a guide for the eye. Here, (3 is obtained 
by averaging the slopes of the individual curves, (b) Log-log 
plot of the power spectral density (PSD) against the angu¬ 
lar frequency (cj) corresponding to dynamical variables of a 
tagged particle (N = 512): (i) the x-component of the veloc¬ 
ity {vx{t)), and (ii) the x-component of position in the COM 
frame (x'(t)). The lengths of the time series’ in both cases are 
equal to 2^® time steps. The straight lines indicate power-law 
fits for the PSD curves. The solid straight line has a neg¬ 
ative slope, = 1.08 and the dashed line has a negative 
slope, ax = 2.21. (c) Illustration of the mechanism through 
which fluctuations in forces occur. Refer to main text for 
details, (d) Lyapunov spectrum showing the AN exponents 
with N = 512. The spectrum is obtained after a convergence 
time lO’^ time steps. The inset plot shows the five largest 
Lyapunov exponets for values of m = 1 (square), 2 (triangle 
up), 3 (circle), 4 (triangle down) and 5 (diamond). 


velocities be aligned and of magnitude vq. This configu¬ 
ration results in A, experiencing attractive forces (shown 
with solid black arrows), small in magnitude and directed 
towards individual particles. The resultant force Fq^a is 
an average vector, and is expected to be small, as well. 
Now, consider a perturbation, displacing the particle D 
by small amount such that its new position is above the 
line La- The configuration would result in E replacing 
D in the set Qa- Due to the location of A, the new Fq,a 
would be of larger magnitude and would have a differ¬ 
ent direction. In fact, the direction of force is expected 
to further increase the separation with D and would not 
allow D to get restored as neighbour. The above mecha¬ 
nism prevents the particles to settle down into a ‘frozen’ 
structure in the COM frame, when beginning from arbi¬ 
trary initial conditions. As to increases, because of the 
presence of a larger number of neighbours, averaging of 
the force is better and the characteristic fluctuation is 
lesser. For to < 3, certain stable structures may be re¬ 
tained when taken as initial conditions. However, such 


















5 


structures are not observed for arbitrary initial condi¬ 
tions. 

The fact that perturbations are able to proliferate in 
the system, generally implies sensitivity to initial con¬ 
ditions. This prompted us to calculate the Lyapunov 
spectrum for the dynamics. We use the Gram-Schmidt 
orthonormalization procedure [23, 24]. The converged 
spectrum for N = 512 and m = 3 with 47V exponents is 
shown in Fig. 5(d). In the inset we compare this with 
the spectra for other values of m by plotting the first 
five Lyapunov exponents. For all the values of m, the 
system has positive exponents and hence indicate expo¬ 
nential instability. For m = 1 and m = 2, the motion of 
particles remain confined to a region (V ~ 0), and the 
largest Lyapunov exponents (LLEs) are 1.72 and 1.42, 
respectively. This suggests that the motion is chaotic. 
However, for larger values of m the motion is unbounded 
and LLEs are smaller (0.28, 0.17, and 0.08 for to = 3, 
4 and 5, respectively). We believe that for to = 3 the 
flock is poised at the edge of chaos, which is generally 
assumed to be the case with SOC [25]. The nature of 
the dynamics is very similar to that of the train model 
[26], where the self-organized critical states were shown 
to have a positive LLE. 

It turns out that the cohesion between particles is a 
key ingredient for the development of long-ranged corre¬ 
lations. The calculation of the force, in case of individ¬ 
ual particles, implicitly depends on the knowledge of the 
COM of the system. This rule allows the flock to be cohe¬ 
sive unit even in the presence of fluctuations and the flock 


is always in the zero-density limit. A simple topological 
rule may not be sufficient to prevent the fragmentation of 
the flock into smaller clusters and subsequent expansion 
in the space, when simulated in open boundary condi¬ 
tions. We would like to argue that the case in real flocks 
may not be completely different. Birds may actually de¬ 
tect the density of individuals over a region and might 
have a preference of joining regions of higher densities, 
and, in effect, getting pulled towards the COM. Bacte¬ 
rial aggregates are studied over a quasistationary state, 
where a finite density is guaranteed. Also, note, that in 
the steady state of the model, the actual forces appear 
locally on scales much smaller than the flock size. 

In conclusion, we provide an SPP model with a mech¬ 
anism through which critical fluctuations can appear in¬ 
side flocks in the presence substantial amount of order 
and in the absence of external noise. Noise in the flock 
is rather internal, where, force fluctuations are caused by 
reshuffling of neighbourhoods and vice versa. We have 
also included the concept of topological distance while 
defining neighbourhood of particles. In future, we plan 
to include the effect of noise and external fields in our 
model, and characterize the response of the flock. Also, 
it would be interesting to construct such a model in the 
evolutionary paradigm [27] where the susceptible state 
naturally evolves as stable strategy. 
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